/*
 * CTDR_interaction.cpp
 *
 *  Created on: 2 déc. 2010
 *      Author: mbenoit
 */

#include "CTDR_Interaction.h"
#include "CTDR_Electron.h"
#include "CTDR_Hole.h"

/*************************************************************************
 * . Class CTDR_Interaction				                 *
 *   Written by : Mathieu Benoit                                         *
 *   m.benoit@umontreal.ca						 *
 *************************************************************************/
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <string>
#include "global.h"
using namespace std;

//////////////////////////////////////////////////////////////////////////
// Source file : class CTDR_Interaction					//
// This class represent a cloud of charge element of a gamma 		//
// CTDR_Interaction in CZT.It contains all the function to create 	//
// and simulate the behavior of an carrier cloud generated by a 	//
// gamma event. The cloud consiste of CTDR_Electron objects,		//
// wich all evolves in time from their generation to their collection. 	//
// 									//
// Constant to set in source before compilation				//
// "N" represent the number of charge element (class CTDR_Electrons)//
// used to reproduce the real charge cloud generated by a gamma 	//
// CTDR_Interaction. The higher the number, he more accurately the 	//
// behavior is simulated. However , higher number of charge element also//
// mean longer computation time. CTDR_Fields and weighting potentials 	//
// must be extracted from solution files obtained with poisson equation //
// solvers for each charge at each time step.				//
// To simulate signals and to focus on one event or one multi-hit, 	//
// a high number of charge is reccomended. But for simulation of a large//
// number of events, for example to obtain scatter plots, a small number//
// of charge works fine and give much more reasonable. As computer 	//
// become faster, these type of simulation will gain much precision as 	//
// they require a lot of brutal force.These simulation can be performed //
// on a single-processor , but the time needed to obtain readable data 	//
// is quite long,(days or weeks). These simulations can easily be 	//
// performed in parallel on many processor.The average parallel		//
// computing facility these days could obtain good results quite fast 	//
// for simulation with a high number of charge elements. This code 	//
// can easily scale to the demand and to the ressource at your 		//
// disposition to give  good results in a maximum of situations.	//
//									//
// zmax is the maximum height of the CZT block 				//
//////////////////////////////////////////////////////////////////////////
//const int N=10;
//ClassImp(CTDR_Interaction)
//______________________________________________________________________________
CTDR_Interaction::CTDR_Interaction(double x, double y, double z,double Energie,double sigma,double t0,int Nelec,int Q){
//Constructor for the CTDR_Interaction , each CTDR_Interaction built represent
//one CTDR_Interaction site with deposited energy "Energie", and the cloud
//begins its simulation with a radius of sigma. It has been found that
//Electron-hole separation process is fast and can be neglected.. ,


	double xtemp; // Centroid position in X
	double ytemp; // Centroid position in Y
	double ztemp;// Centroid position in Z
	energie=Energie;  // Interaction's energy
	this->Q=Q;
	N=Nelec; //The number of charge elements used for the simulation
	t=t0;
	Dt=0.001; // Dt in nanoseconds

	// The charge elements are distibuted according to a gaussian of sigma in each dimension
	// Around the CTDR_Interaction site position

	cout << TString::Format("[CTDR_Interaction] Generating interaction point with energy %f at (x,y,z)=%f,%f,%f",Energie,x,y,z) << endl;

	for(int i=0;i<N;i++){

		xtemp=gCTDR.rand->Gaus(x,sigma);
		ytemp=gCTDR.rand->Gaus(y,sigma);
		ztemp=gCTDR.rand->Gaus(z,sigma);

		//We set the mobility, 0.1 for CTDR_Electron in CZT, and the amount of real charge
		// attributed to each charge elements

		if(Q<0){
			CTDR_Electron carrier;
			carrier.SetValeurInitial(xtemp,ytemp,ztemp,Default_conv_Electron_Mobility,-1);
			carrier.SetD(Default_conv_Electron_D);
			charges.push_back(carrier);

		}
		else{
			CTDR_Hole carrier;
			carrier.SetValeurInitial(xtemp,ytemp,ztemp,Default_conv_Hole_Mobility,1);
			carrier.SetD(Default_conv_Hole_D);
			charges.push_back(carrier);


		};
		};
};


CTDR_Interaction::CTDR_Interaction(double x, double y,double Energie,double sigma,double t0,int Nelec,int Q){
//Constructor for the CTDR_Interaction , each CTDR_Interaction built represent
//one CTDR_Interaction site with deposited energy "Energie", and the cloud
//begins its simulation with a radius of sigma. It has been found that
//Electron-hole separation process is fast and can be neglected.. ,


	double xtemp; // Centroid position in X
	double ytemp; // Centroid position in Y
	double ztemp;// Centroid position in Z
	energie=Energie;  // Interaction's energy
	this->Q=Q;
	N=Nelec; //The number of charge elements used for the simulation
	t=t0;
	Dt=0.001; // Dt in nanoseconds

	// The charge elements are distibuted according to a gaussian of sigma in each dimension
	// Around the CTDR_Interaction site position

	//cout << TString::Format("[CTDR_Interaction] Generating interaction point with energy %f at (x,y,z)=%f,%f,%f",Energie,x,y,z) << endl;

	for(int i=0;i<N;i++){
//		Rtemp=gCTDR.rand->Uniform(0,1);
//        R2temp=gCTDR.rand->Uniform(0,1);
//        xtemp=x+sigma*sin(2*pi*Rtemp)*sqrt(-2*log(R2temp));
		//xtemp=gCTDR.rand->Gaus(x,sigma);
		  xtemp=x;

//        Rtemp=gCTDR.rand->Uniform(0,1);
//        R2temp=gCTDR.rand->Uniform(0,1);
//        ytemp=y+sigma*sin(2*pi*Rtemp)*sqrt(-2*log(R2temp));
		//ytemp=gCTDR.rand->Gaus(y,sigma);
		  ytemp=y;
//        Rtemp=gCTDR.rand->Uniform(0,1);
//        R2temp=gCTDR.rand->Uniform(0,1);
//		ztemp=z+sigma*sin(2*pi*Rtemp)*sqrt(-2*log(R2temp));
		//ztemp=gCTDR.rand->Gaus(z,sigma);
		ztemp=gCTDR.rand->Uniform(0,zmax);
		//We set the mobility, 0.1 for CTDR_Electron in CZT, and the amount of real charge
		// attributed to each charge elements
		if(Q<0){
			CTDR_Electron carrier;
			carrier.SetValeurInitial(xtemp,ytemp,ztemp,Default_conv_Electron_Mobility,-1);
			carrier.SetD(Default_conv_Electron_D);
			charges.push_back(carrier);

		}
		else{
			CTDR_Hole carrier;
			carrier.SetValeurInitial(xtemp,ytemp,ztemp,Default_conv_Hole_Mobility,1);
			carrier.SetD(Default_conv_Hole_D);
			charges.push_back(carrier);


		};
		};
};
//______________________________________________________________________________
CTDR_Interaction::CTDR_Interaction(vector<double> x, vector<double>  y, vector<double> z,double Energie,double sigma,double t0, int Q)
{
//Constructor for the CTDR_Interaction , each CTDR_Interaction built represent
//one CTDR_Interaction site with deposited energy "Energie", and the cloud
//begins its simulation with a radius of sigma. It has been found that
//Electron-hole separation process is fast and can be neglected.. ,

	energie=Energie;  // Interaction's energy
	N=x.size(); //The number of charge elements used for the simulation
	t=t0; // t0 in nanoseconds
	Dt=1.0; // timestep in nanosecond
	//double dE= Energie/N;
	// The charge elements are distibuted according to a gaussian of sigma in each dimension
	// Around the CTDR_Interaction site position

	cout << TString::Format("[CTDR_Interaction] Generating custom interaction point") << endl;

	for(int i=0;i<N;i++){

		CTDR_Carrier carrier;
		if(Q<0)
		{
			CTDR_Electron carrier;
			carrier.SetValeurInitial(x.at(i),y.at(i),z.at(i),Default_conv_Electron_Mobility,Q);
			carrier.SetD(Default_conv_Electron_D);
			charges.push_back(carrier);

}
		else{
			CTDR_Hole carrier;
			carrier.SetValeurInitial(x.at(i),y.at(i),z.at(i),Default_conv_Hole_Mobility,Q);
			carrier.SetD(Default_conv_Hole_D);
			charges.push_back(carrier);

};
		};
};


CTDR_Interaction::CTDR_Interaction(double track[6], double Energie,int Nelec,double Q)
{

	double xtemp,ytemp,ztemp;
	energie=Energie;  // Interaction's energy
	this->Q=Q;
	N=Nelec; //The number of charge elements used for the simulation
	t=0;
	Dt=0.001; // Dt in nanoseconds

	double dx= (track[1]-track[0]);
	double dy= (track[3]-track[2]);
	double dz= (track[5]-track[4]);
	double dydz=dy/dz;
	double dxdz=dx/dz;

	cout << TString::Format("[CTDR_Interaction] Generating track with entry point at (x,y,z)=%f,%f,%f and exit point at (x,y,z)=%f,%f,%f",track[0],track[2],track[4],track[1],track[3],track[5]) << endl;


	for(int i=0;i<N;i++){

		ztemp=gCTDR.rand->Uniform(track[4],track[5]);
		xtemp=track[0]+ (ztemp-track[4])*dxdz;
		ytemp=track[2]+ (ztemp-track[4])*dydz;



		if(Q<0)
		{
			CTDR_Electron carrier;
			carrier.SetValeurInitial(xtemp,ytemp,ztemp,Default_conv_Electron_Mobility,-1);
			carrier.SetD(Default_conv_Electron_D);
			charges.push_back(carrier);
		}
		else
		{
			CTDR_Hole carrier;
			carrier.SetValeurInitial(xtemp,ytemp,ztemp,Default_conv_Hole_Mobility,1);
			carrier.SetD(Default_conv_Hole_D);
			charges.push_back(carrier);
		}
	};

	//cout << "[CTDR_Interaction] il y a " << charges.size() << endl;

};



//______________________________________________________________________________
bool CTDR_Interaction::iteration(CTDR_Field &A)
{
//This function perform a time step for all the charge element in the cloud

       bool I=0;
       this->ComputeRms();
	   #pragma omp parallel for
       for(int i=0;i<N;i++){
    	double z_border= charges.at(i).GetZ();
    	if(z_border<zmax && z_border>0){
    	//If charge have not reached the anode
		//Charges cannot go lower than the cathode
		//We do a time step for each charge element
    	A.SetChamp(charges.at(i).GetX(),charges.at(i).GetY(),charges.at(i).GetZ());

    	if(charges.at(i).GetQ()<0)charges.at(i).MobilityElectron(A.GetE());
    	else charges.at(i).MobilityHole(A.GetE());

    	charges.at(i).Iterateur(Dt,int(energie/4.4),rmsX,rmsY,rmsZ,A); };

    	if(z_border<=0){charges.at(i).Collecte(0);};

    	if(z_border>zmax){charges.at(i).Collecte(zmax);};


	//If no charge are free, I is 1, else its 0
        if(charges.at(i).GetEcollecte()==Libre){I=1;};
		};

		t+=Dt; //time advance

		//If all charge are collected , we return the signal to stop for this CTDR_Interaction to the global simulation
		if(I==0){return 1;}
		else {return 0;};
		//cout << t << endl;

};

//______________________________________________________________________________
void CTDR_Interaction::ComputeRms()
{
//This function computes the RMS in x,y,z of the charge cloud, the centroid is also calculated
//These data are passed at each time step the the charge elements to set-up their random walk

	double carrex=0,carrey=0,carrez=0,moyx=0,moyy=0,moyz=0;
	ErreurMoy=0;
	int nlibres=0;
	for(int i=0;i<N;i++) {
		carrex+=pow(charges.at(i).GetX(),2);
		carrey+=pow(charges.at(i).GetY(),2);
		carrez+=pow(charges.at(i).GetZ(),2);
		moyx+=charges.at(i).GetX();
		moyy+=charges.at(i).GetY();
		moyz+=charges.at(i).GetZ();
		if(charges.at(i).GetEtat()==libre){
				ErreurMoy+=charges.at(i).GetErreur();
				nlibres++;
				};
		};
	rmsX=sqrt(fabs(carrex/N - pow(moyx/N,2)));
	rmsY=sqrt(fabs(carrey/N - pow(moyy/N,2)));
	rmsZ=sqrt(fabs(carrez/N - pow(moyz/N,2)));
	MoyZ=moyz;

	ErreurMoy/=nlibres;
	SetDt();
	//cout << "rms x,y,z: " << rmsX << ' ' << rmsY << ' ' << rmsZ << " t=" << t << endl ;



	cout << TString::Format(" t=%03.3f (avg X,avg Y,avgZ)=(%05.2f,%05.2f,%05.2f) ",t,moyx/N,moyy/N,moyz/N) ;

};

//______________________________________________________________________________
double CTDR_Interaction::GetRmsX()
{
//Access method to RMS

	return rmsX;
};

//______________________________________________________________________________
double CTDR_Interaction::GetRmsY()
{
//Access method to RMS
	return rmsY;
};

//______________________________________________________________________________
double CTDR_Interaction::GetRmsZ()
{
//Acces method to RMS
	return rmsZ;
};

//______________________________________________________________________________
double CTDR_Interaction::GetMoyZ()
{
//Acces method to centroid position in Z
	return MoyZ/N;
};

//______________________________________________________________________________
double CTDR_Interaction::GetQx()
{
//Acces method to centroid position in Z

	return Q0x;
};

//______________________________________________________________________________
double CTDR_Interaction::GetQy()
{
//Acces method to centroid position in Z

	return Q0y;
};

//______________________________________________________________________________
void CTDR_Interaction::SetQx(double Q)
{
//Modifier function of the initial weight of the charge cloud. Data used in signal computation

	Q0x=Q;
};

//______________________________________________________________________________
void CTDR_Interaction::SetQy(double Q)
{
//Modifier function of the initial weight of the charge cloud. Data used in signal computation

	Q0y=Q;
};

//______________________________________________________________________________
void CTDR_Interaction::SetDt()
{

	//cout << ErreurMoy <<endl;

	if(isnan(ErreurMoy)){Dt=tup;}
	else if(ErreurMoy > Target){ Dt*=0.9;}
	else if(ErreurMoy < Target){ Dt*=1.1;};


	if(Dt<tlow) Dt=tlow;
	if(Dt>tup)  Dt=tup;

};

//______________________________________________________________________________
int CTDR_Interaction::GetN()
{
// Access method for N

	return N;
};

//______________________________________________________________________________
CTDR_Carrier CTDR_Interaction::GetCharge(int i)
{
// Access method for Electrons object
	return charges.at(i);
};

//______________________________________________________________________________
double CTDR_Interaction::Gett()
{
// Access method for time
	return t;
};

//______________________________________________________________________________
void CTDR_Interaction::WritePositions(ofstream &sortie)
{
// This function write the position of the charge elements to a file


	for(int i=0;i<N;i++){
		sortie << charges.at(i).GetX() << ' ' << charges.at(i).GetY() << ' '
		<< charges.at(i).GetZ() << ' ' << charges.at(i).GetQ() << endl;
		};
};
